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SUMMARY 


Unsustainable exploitation of natural resources is 
increasingly affecting the highly biodiverse tropics 
[i, 2]. Although rapid developments in remote 
sensing technology have permitted more precise 
estimates of land-cover change over large spatial 
scales [3-5], our knowledge about the effects of 
these changes on wildlife is much more sparse 
[6, 7]. Here we use field survey data, predictive den- 
sity distribution modeling, and remote sensing to 
investigate the impact of resource use and land-use 
changes on the density distribution of Bornean 
orangutans (Pongo pygmaeus). Our models indicate 
that between 1999 and 2015, half of the orangutan 
population was affected by logging, deforestation, 
or industrialized plantations. Although land clear- 
ance caused the most dramatic rates of decline, it 
accounted for only a small proportion of the total 
loss. A much larger number of orangutans were lost 
in selectively logged and primary forests, where rates 
of decline were less precipitous, but where far more 
orangutans are found. This suggests that further 
drivers, independent of land-use change, contribute 


updates 
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to orangutan loss. This finding is consistent with 
studies reporting hunting as a major cause in orang- 
utan decline [8-10]. Our predictions of orangutan 
abundance loss across Borneo suggest that the 
population decreased by more than 100,000 individ- 
uals, corroborating recent estimates of decline [11]. 
Practical solutions to prevent future orangutan 
decline can only be realized by addressing its com- 
plex causes in a holistic manner across political 
and societal sectors, such as in land-use planning, 
resource exploitation, infrastructure development, 
and education, and by increasing long-term sustain- 
ability [12]. 


RESULTS 


Bornean Orangutan Field Survey Data 

To model Bornean orangutan (Pongo pygmaeus) density distri- 
bution and derive metapopulation abundances, we compiled 
orangutan field surveys. Estimates of orangutan density and 
abundance are usually derived from the observation of their 
nests [13, 14] on line transects [15]. A total of 36,555 orangutan 
nests were observed on 1,491 ground and 252 aerial transects 
that were surveyed between 1999 and 2015 throughout the 
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Bornean orangutan range, with a total survey effort of 4,316 km 
(ground: 1,388 km; aerial: 2,928 km) and a median of 86 tran- 
sects (interquartile range [IQR]: 28-156 transects) per year. The 
cumulative area of land surveyed contained 1,234 km?. During 
the study period, the average yearly encounter rate significantly 
decreased from 22.5 to 10.1 nests/km (parameter estimate = 
—0.06, SE = 0.02, Z = —2.25, p = 0.04; the model contained 
the log-transformed mean nest encounter rate per year as 
response, weighted by the number of transects per year and 
the year as predictor). 


Estimating Change in Bornean Orangutan Density 
Distribution 

We built a predictive density distribution model to estimate 
Bornean orangutan abundance. The full model included survey 
year, climate, habitat cover, and human threat predictor 
variables (see STAR Methods and Key Resources Table) and 
explained orangutan density significantly better than the null 
model including only the intercept (likelihood ratio test, 
X? = 1,440, df = 13, p « 0.001). Mean temperature and lowland 
and peatswamp forest cover had a significant positive relation- 
ship with orangutan density (Figure S1; Table S1). Study year, 
rainfall variability, montane forest cover, and human population 
density negatively affected orangutan density (Figure S1; 
Table S1). Intermediate levels of rainfall in dry months were 
related to higher densities of orangutans. The cover of defor- 
ested areas around transects was slightly positively correlated, 
but its confidence limits included zero. Topsoil organic carbon 
content, estimate of orangutan killing, and percentage of the 
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population with hunting taboos were not significantly correlated 
with orangutan density. 

With the aim of minimizing model uncertainty in spatial model 
predictions, we used multi-model inference and evaluated all 
possible combinations of covariates included in the full model 
(Table S1). The complete set of all fitted models was then 
used to estimate the orangutan density distribution across the 
range. The estimated distribution was mapped to metapopula- 
tions delineated by experts at the Population and Habitat 
Viability Assessment Workshop (PHVA) for Bornean orangu- 
tans. In this context, the term "metapopulation" was used to 
identify larger entities that are bound by dispersal barriers, 
such as rivers, major roads, and areas without forests, and 
that include one or more orangutan subpopulations. Only 38 
out of 64 identified metapopulations retained more than 100 
individuals and can thus be considered to contain viable sub- 
populations [16]. 

The three largest metapopulations were found in Kalimantan, 
the Indonesian part of Borneo, and have experienced a strong 
decline over the studied 16-year period (Figure 1). 

Western Schwaner, the largest metapopulation, lost an esti- 
mated 42,700 individuals (9596 confidence interval [CI]: 
12,700-73,400) since 1999, with 40,700 (95% Cl: 30,000- 
57,200) remaining in 2015. The second-largest population, 
Eastern Schwaner, lost 20,100 individuals (95% Cl: 7,200- 
33,500) and was estimated to contain 16,800 (9596 CI: 12,100- 
23,100) in 2015. In Karangan, the third-largest population, 
8,200 individuals (9596 CI: 1,900-15,400) were lost and 9,000 
(5,900-14,200) remained in 2015. The total estimated loss of 


1 - Western Schwaner 


- 
O a 
o o 


o o 
o o 
o o 


o 
eo 
e 
e 
e 


Orangutan abundance 


2000 2005 2010 


year 


2 - Eastern Schwaner 


Orangutan abundance 


2000 


2005 2010 


year 








2015 2020 2050 





2015 2020 2050 


3 - Karangan 


2000 2005 2010 


year 


2015 2020 2050 


Figure 1. Abundance of the Three Largest Orangutan Metapopulations between 1999 and 2015 and Projected Abundance for 2020 and 2050 
Orangutan abundance was estimated for the three largest metapopulations with a multi-model approach over the study period (1999 to 2015). Estimates of future 
orangutan abundance were based on forest cover projections for 2020 and 2050 by Struebig et al. [17] and are indicated by a dashed line. Shaded areas and error bars 
represent the 95% confidence intervals. On the y axes, the number “10,000” is highlighted in blue to show the scale difference between the three populations. The map 
shows all identified metapopulations in gray. The three largest metapopulations are indicated by their color. State labels are as follows: Br, Brunei; Sb, Sabah; and Sk, 
Sarawak in Malaysia; WK, West; EK, East; NK, North; SK South; and CK, Central Kalimantan in Indonesia. See also Figures S1 and S2 and Tables S1-S3. 


Bornean orangutans between 1999 and 2015 amounted to 
148,500 individuals (95% Cl: 48,100-252,300). 

We used predictions of forest cover from Struebig et al. [17] for 
2020 and 2050 to project future orangutan decline (Figures 1 
and 2). To this end, we assumed that orangutans cannot survive 
in areas without tree cover. The orangutan abundance in the 
three largest populations was projected to drop further and 
reach 31,100 individuals (9596 Cl: 22,500-44,000) in the Western 
Schwaner metapopulation area, 14,700 individuals (9596 CI: 
9,600-19,600) in Eastern Schwaner, and 6,100 individuals 
(95% Cl: 3,800-10,000) in Karangan by 2050. The total future 
loss for all metapopulations was projected to be 45,300 (95% 
Cl: 33,300-63,500). This projected future decline is only based 
on the direct consequence of habitat loss. It does not consider 
the effects of orangutan killing for food and in conflict and is 
therefore most likely an underestimate. All estimates are 
rounded to the nearest hundred. 


Linking Remotely Sensed Resource Use and Density 
Distribution 
To identify possible causes for the estimated orangutan loss, we 
compared absolute abundance and density from the beginning 
and the end of the survey period between land-use types and as- 
sessed differences in change over time. We differentiated areas 
in which resource use had altered the environment and areas in 
which land-use remained unaltered during the study period. For 
land-use changes, we considered deforestation, conversion to 
industrial plantations (oil palm and paper pulp), and selective 
logging in natural forests. As stable land-use, we considered pri- 
mary and montane primary forest, regrowth forests, industrial 
plantations established prior to the study period, and “other,” 
comprising non-forest areas. 

By 2015, 50% of the orangutans estimated to have occurred 
on Borneo in 1999 were found in areas in which resource use 
had altered the environment. A comparison of distinct regions 
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Figure 2. Spatial Distribution of Estimated Orangutan Densities on Borneo for the Year 1999 and 2015 and Projections to 2020 and 2050 
Bornean orangutan density per 1 km? in the beginning and the end of the study period and for 2020 and 2050. Between 1999 and 2015, high-density areas (dark 
green) disappeared, and medium-density areas (light green) declined. Low-density areas (beige and purple) expanded. Future estimates are based on projected 
forest loss [17]; therefore, map representations between model estimates and future projections differ. Areas in which forest was projected to be lost also lose the 
resident orangutans. Hence, maps between 2015 and 2020 seem to lose many fragments inhabited by orangutans, but they already had low density before. 
Between 2020 and 2050, further areas were projected to lose forest, but the loss is less visible. See also Figures S1 and S2 and Tables 81-53. 


revealed that 5096, 6096, and 1096 of the orangutans were 
affected by transformation into industrial oil palm or paper pulp 
plantations, deforestation, or selective logging in Kalimantan, 
Sabah, and Sarawak, respectively. Rates of orangutan decline 
were highest in areas deforested or converted to plantations 
(6396-7596 loss) in both Kalimantan and Sabah (Figure 3). In Sar- 
awak, there were almost no industrial plantations and deforested 
areas within the orangutan metapopulation range, together 
affecting only 0.496 of the area and 296 of the orangutan popula- 
tion. Industrial plantations and deforestation contributed 796 
(Kalimantan), 296 (Sabah), and less than 196 (Sarawak) to the 
overall estimated loss of orangutans in each of the three regions. 

Both Kalimantan and Sabah had the highest orangutan abun- 
dance in selectively logged forests, followed by primary forest. In 
Sarawak, the highest orangutan abundance was found in pri- 
mary forests. The rate of orangutan decline across the three re- 
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gions and these two land-use classes was less precipitous, but 
still high (4996-5696). The loss of orangutans in primary and 
selectively logged forests between 1999 and 2015 accounted 
for 6796 of the total loss in Kalimantan (93,000 individuals; 
9596 Cl: 26,500-162,300), 7296 in Sabah (6,100 individuals; 
9596 Cl: 2,400-10,000), and 83% of the total loss in Sarawak 
(900 individuals; 95% Cl: 250-1,600). 


DISCUSSION 


The unsustainable use of natural resources has caused a 
dramatic decline of Bornean orangutans. Only 38 out of 64 
remaining metapopulations have more than 100 individuals, 
the assumed threshold for viability of Bornean orangutan 
populations [16]. Our findings suggest that more than 100,000 in- 
dividuals have been lost in the 16 years between 1999 and 2015. 
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Figure 3. Linking Remotely Sensed Resource Use and Density Distribution 

Percent area affected by resource use in orangutan metapopulations during the study period, forest and non-forest classes (pie charts), their spatial distribution 
(map), and total orangutan abundance and its change between the first study year (1999) and last study year (2015) (bar charts). Total metapopulation areas per 
province in square kilometers are given in the lower-right corner of the pie charts. Areas had been transformed into plantations (oil palm and paper pulp), de- 
forested, or selectively logged between 1999 and 2015; were covered with forest (regrowth, primary or montane primary forest); were plantations already before 
the study period; or were another unspecified non-forest class. The percent orangutan abundance loss in comparison to 1999 is highlighted in rectangles in the 
bar charts. The “*” indicates the absence of orangutans in the respective category. The error bars indicate the 95% confidence interval. On the x axes, the number 
“2000” is highlighted in blue to show the scale differences between the three areas. See also Figure S3. 
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All three analytical approaches employed in this study, based on 
field survey data, spatial covariate modeling, and remote 
sensing, corroborated the concluded impact of resource use 
and resulting decline of Bornean orangutans. The results are 
also very consistent with the genetic signature of a recent 
collapse found in an orangutan population in Sabah [18] and 
evidence of large annual losses of orangutans through hunting 
and conflict killing in Kalimantan [8-10]. Our results substantiate 
the percentage loss estimated by Santika et al. [1 1] and reinforce 
the recent uplisting of the Bornean orangutan as Critically 
Endangered on the IUCN Red List [19]. The numbers reported 
here are larger than past estimates [11], but they are in line 
with findings reported for other great ape taxa [20-23]. 

We have established the density distribution of Bornean 
orangutans with a model-based approach that uses the relation- 
ships between predictor variables and observed orangutan 
abundance to predict abundance for unsurveyed sites. These 
predictions are useful for deducing trends at the regional to 
landscape scale [24], but they may be limited at a local scale, 
where additional demographic and behavioral drivers can influ- 
ence orangutan density distribution, e.g., ranging behavior in 
response to local food resources or conspecifics. Thus, our find- 
ings reveal patterns at large spatial scales, but great care should 
be taken when inferring from predictions at specific sites. 

Another aspect of our study that requires critical assessment 
is the inference of orangutan abundance from nest counts. 
Nest decay time, an essential parameter to translate nest density 
into orangutan density, varies between survey sites. Although 
factors like rainfall, wood density, and complexity of nest 
architecture are known to influence nest decay time [13, 25, 
26], additional variability in decay time between sites is not fully 
understood [27]. We addressed this issue by using all available 
datasets on orangutan nest decay, comprising information on 
the lifespan of more than 1,000 nests (see STAR Methods) 
across Borneo. If our findings of orangutan decline were an arti- 
fact of severely biased nest decay times, this would require nest 
decay time to have halved over the course of the study period. 
However, we found no indication of this and so do not consider 
this to be a limitation of our study. 

Contrary to our expectations, the model coefficient for defores- 
tation indicated a slightly positive relationship between deforesta- 
tion in years prior to the survey and orangutan abundance. There 
are several possible explanations for this observation, suggesting 
that the model coefficient does not capture a causal relationship. 
First, surveys tend to be biased toward areas with known orang- 
utan occurrence. Thus, our dataset possibly lacks sufficient vari- 
ance for detecting the true impact of deforestation on orangutan 
density. Second, some studies have suggested that the number 
of orangutans in areas adjacent to deforested areas is temporally 
inflated, due to the displacement of individuals and subsequent 
refugee crowding [28, 29]. Third, high dietary flexibility allows 
orangutans to be resilient in the face of some levels of disturbance 
[30, 31]. This may delay the effects of deforestation on the 
observed density for several years, before populations eventually 
start to decline [28]. Irrespective of this, when we compare spatial 
model predictions and remotely sensed land-use change, the 
highest rates of orangutan decline were detected in areas with 
habitat removal (deforestation and conversion to industrial planta- 
tions). This shows that the predictive density distribution model 
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has indirectly captured the deleterious effects of deforestation 
on orangutan abundance. Our finding suggests that deforestation 
and industrial oil palm and paper pulp plantations are responsible 
for about 996 (14,000 individuals) of the total loss of orangutan 
abundance. Whereas in the early years of the study it was mainly 
degraded land with low orangutan density that was converted to 
industrial plantations, after 2005 the conversion of forests to oil 
palm plantations has been increasing dramatically [32]. Some 
studies have suggested that orangutans can occur in oil palm or 
paper pulp plantations, when these are managed well and adja- 
cent forest fragments are maintained [33-35]. However, it is 
unclear whether this is just a transient effect or whether orangu- 
tans can indeed persist over the long term [33-35]. 

The highest orangutan abundances were found in selectively 
logged forests in Kalimantan and Sabah and in primary forests 
in Sarawak. This finding is consistent with studies reporting 
that orangutans can occur in selectively logged or regenerating 
logging concessions, depending on the type and intensity of 
logging operations [36-39]. Consequently, successful orangutan 
conservation is necessarily situated in multi-functional land- 
scapes [36, 40] and recognizes the importance of degraded 
and logged forests as well as forest fragments in plantation 
matrices [33, 34]. 

Effective partnerships with logging companies, whose conces- 
sions harbor the majority of orangutans, are essential to curb 
orangutan loss [41]. Similarly, partnerships with oil palm and pa- 
per pulp producers are important to promote best practice guide- 
lines for management [33, 35, 42]. Such partnerships have already 
been reported, e.g., by Meijaard et al. [43], and could potentially 
provide co-benefits for biodiversity conservation in general [37]. 
The Roundtable on Sustainable Palm Oil (RSPO) and the Forest 
Stewardship Council (FSC) are examples of certification schemes 
that incentivize these partnerships, by enabling consumers to 
favor responsible natural resource management [42]. 

The pervasive decline of orangutans in more intact habitat is 
consistent with various studies identifying hunting as the main 
driver of biodiversity loss inthe tropics [44, 45], including Southeast 
Asia [2]. More specifically, our observation is supported by the re- 
sults of extensive interview surveys in Kalimantan that show that, 
per year, on average 2,256 orangutans were hunted or killed due 
to conflict with humans [8-10]. The estimate of orangutan killing 
in the model is based on a Borneo-wide projection of hunting pres- 
sure derived from these interview surveys [10]. In the model, this 
predictor did not show an influence on orangutan density. It is 
possible that our survey dataset lacks sufficient variance for de- 
tecting the impact of killing on orangutan density or the available 
predictor layer does not well represent the actual pressure, as 
the relationship between density and hunting or killing is very com- 
plex and non-linear. Human population density, on the other hand, 
had a significant negative influence on orangutan densities in the 
model and may have already captured the effect of orangutan 
killing. Orangutans are also present in the national and international 
wildlife trade. Traded orangutans are usually young orphans, and 
for each orphan, adult individuals have been killed [46]. Due to 
thelow reproductive rate ofthe species, even very low offtake rates 
of reproductive females (~1% per year) will drive populations to 
extinction [16, 47]. In the absence of plausible alternative explana- 
tions for the observed loss of orangutans in seemingly intact 
habitats, such as the occurrence of widespread and highly lethal 


infectious diseases as observed among African apes [48], killing is 
the most likely explanation. From this perspective, our prediction 
ofa further loss of 45,300 orangutans over the next 35 years, based 
solely on projections of forest cover change, is most likely an un- 
derestimate. Furthermore, many individuals currently occur in 
fragmented, small populations that are assumed not to be viable 
and will most likely disappear in the near future. 

Knowledge about the density distribution of key species is 
essential to explore the consequences of land-use change, 
exploitation of natural resources, development of infrastructure, 
and climate change. It is also needed to evaluate which conser- 
vation interventions are most effective in reducing decline and 
loss of biodiversity. 

In essence, natural resources are being exploited at unsus- 
tainably high rates across tropical ecosystems, including 
Borneo. As a consequence, more than 100,000 Bornean orang- 
utans vanished between 1999 and 2015. The major causes are 
habitat degradation and loss in response to local to global de- 
mand for natural resources, including timber and agricultural 
products, but very likely also direct killing. Our findings are 
alarming. To prevent further decline and continued local extinc- 
tions of orangutans, humanity must act now: biodiversity conser- 
vation needs to permeate into all political and societal sectors 
and must become a guiding principle in the public discourse 
and in political decision-making processes. 
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http://www.worldclim.org/bioclim; BIO17 


http://www.fao.org/soils-portal/soil-survey/soil-maps- 
and-databases/harmonized-world-soil-database-v12/en/ 


https://ormt-crisp.nus.edu.sg/ 
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global-forest/download_v1.3.html 
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http://www.r-project.org; RRID:SCR_001905 


https://stat.ethz.ch/R-manual/R-devel/library/MASS/ 
html/glm.nb.html 


https://cran.r-project.org/web/packages/car/index.html 
http://www.gdal.org/; RRID:SCR 014396 
https://www.qgis.org/ 

https://www.python.org/; RRID:SCR_008394 


N/A 


Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Maria Voigt (maria. 


voigt@idiv.de). 
METHOD DETAILS 


Study area and orangutan data 


For this study we compiled three types of data: 1) line transect nest count data; 2) nest decay time data; and 3) polygons representing 
areas inhabited by orangutan metapopulations. Bornean orangutan (Pongo pygmaeus) nest count line transect data were compiled 
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from surveys undertaken across Borneo between 1999 to 2015. Researchers reported the number of orangutan nests observed 
along line transects, which were either walked or flown with a helicopter (aerial and ground transects), respectively. The datasets 
were converted to a standard format to include the number of observed nests, total transect length, year of survey, and start and/or 
end coordinates of surveyed transect line. All ground transects with perpendicular distances (ppd) to nests were used for the 
Distance analysis [62] (number of nests = 15,858, 64% of total), to estimate truncation distance and effective strip width (ESW), 
that is, the perpendicular distance from the transect, below which an equal number of nests was missed as seen beyond [14]. For 
the predictive density distribution model we also considered aerial and ground transects without ppd and assumed estimated 
ESW to be representative. The cumulative area of land surveyed was calculated as the transect length multiplied by two times 
the effective strip width, excluding repeat sampling. 

There were only few transects from areas on Borneo in which orangutans are known to be absent. Thus, we added ‘virtual’ 
transects with zero nests randomly to expert-delineated areas of orangutan absence [49] to balance this bias in sampling. For 
each survey year, we set the number of transects in the area of known absences to 50% of the number of surveyed transects in 
the orangutan range in the given year. We tested the effect of varying the number of absence transects (30%, 50% and 80% density 
of surveyed transect), but the model proved to be robust and the resulting orangutan abundance estimate did not differ substantially 
(30% absence density in comparison to 50%: correlation coefficient > 0.99, maximum percent difference = 5.6%; 80% absence den- 
sity in comparison to 50%: correlation coefficient > 0.99, maximum percent difference = 3%; n = 16 years). 

We compiled nest decay information from four sites. For two locations (Sabangau in Central Kalimantan and Lesan in East 
Kalimantan) we used nest decay datasets including information from repeated visits about nest status from construction to disap- 
pearance. The dataset from Lesan included 88 nests, which were visited between February 2005 and September 2006. In Sabangau 
423 nests were visited between July 2001 and April 2011. For two other sites (Kinabatangan, Sabah and Gunung Palung, West Ka- 
limantan) we used information about nest decay time, estimated by Ancrenaz et al. [25] and by Johnson et al. [63]. 

At the PHVA for Bornean orangutans held between the 24" and 27" of May 2016 in Bogor, Indonesia, 41 orangutan experts map- 
ped 64 Bornean orangutan metapopulations [16]. The resulting metapopulation polygons covered areas between 6 and 58, 157 km?, 
amounting to a total area of 333,250 km?. Predictions were extrapolated to this area, and although only a small proportion was actu- 
ally sampled (0.37%), the surveys were distributed well across the area. Only 23% of the metapopulation area was located outside 
the 95% minimum convex polygon of transect locations. 


Predictor variables of orangutan abundance 

We selected predictor variables based on their presumed importance for orangutan ecology, while guaranteeing data availability for 
the whole range and minimizing the correlation between them (Table S2) [24]. The final predictor variable set comprised layers de- 
picting climate (mean daily temperature, yearly variation in rainfall, rainfall in dry months (May - September), habitat (topsoil organic 
carbon content, peatswamp, lowland and lower montane forest cover), and anthropogenic pressures on orangutans (deforestation, 
human population density, orangutan killing estimates, and percent population with religious hunting taboos). The predictor for 
orangutan killing estimates was based on a Borneo wide model of orangutans killed in years prior to interview surveys [8] by Abram 
et al. [10]. We included percent Muslim population as a proxy for the proportion of the population that has hunting taboos, because it 
had been shown that hunting pressure on primates is lower in areas inhabited by a majority of Muslims [9, 64]. 

Before extraction, we reprojected all predictor layers to the Asia South Albers Equal Area Conic, to allow for accurate represen- 
tation of metric distances. The layers were resampled to the same extent, origin and a resolution of 1 km, the coarsest available. 
Nearest neighbor resampling was used for categorical predictors. 

We extracted climate and habitat variables within a radius of 1 km around each transect, resulting in an area of at least 3.14 km?, 
depending on the transect length. This approximates the size of the home range of female orangutans on Borneo and ensures that 
climatic and ecological predictors that have an effect on the population are appropriately represented. Variables indicating 
anthropogenic pressures were obtained within a distance of 10 km, approximating the distance over which human influence is 
most likely (E.M., unpublished data). 

Information about habitat cover was available for three time points (2000, 2010 and 2015 [52, 53]). We used the habitat cover 
information from 2000 for all transects surveyed between 1999 - 2005, the layer from 2010 for all transects surveyed between 
2006 and 2012, and the layer from 2015 for transects sampled in 2013 to 2015. At the time of the analysis, deforestation maps 
were available for each year between 2000 and 2014 [3]. For each transect, we considered the percent area deforested in the years 
prior to the survey in a 10 km-buffer around the transect. 

When the start or the end-point of a transect was unknown, we extracted the predictor variables with a radius of half the transect 
length [sensu 55]. We determined the proportion of each class within the neighborhood for categorical and the mean value for contin- 
uous predictor variables. 

We repeated the extraction for a 1 x 1 km grid covering the metapopulation areas, to enable the estimation of orangutan 
abundance over the whole range. It was visually verified that all predictors had an approximately symmetrical distribution, and human 
population density was subsequently log-transformed. We also ensured that the range of variable values extracted for the transect 
observations was broad enough to meaningfully allow prediction to the range of values extracted for the metapopulation areas by 
comparing the distribution of both. We found that the majority of predictors covered more than 7596 of the predictor space to which 
estimates were extrapolated. The exceptions were the predictors deforestation (6396 cover of sampled predictor range), mean 
temperature (5096 cover) and human population density (» 196 cover). For the predictor mean temperature the low values were 
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not included. These occur in high elevation areas, which were sampled less as they are difficult to access and harbor fewer orang- 
utans [28]. The surveys also did not include areas with high human population density. As the density of orangutans decreases to zero 
in high elevation areas and areas with high human population density, the extrapolation error cannot become large. Thus, we did not 
consider the low coverage for these predictors to be a limitation. The cover of predictor values was at most 3% lower, when excluding 
the absence transects, except for rainfall variability. For this predictor, the absence transects increased the cover of predictor values 
by 19%. Finally, all predictors were standardized to a mean of zero and a standard deviation of one to facilitate the comparison of 
model parameters [65]. 


Future orangutan abundance 

We used information about remaining forest cover on Borneo projected for 2020 and 2050 from Struebig et al. [17, 41] together with 
the orangutan density distribution estimated for 2015 and predicted orangutan distribution 5 and 35 years after the last study year. 
Assuming that orangutans will not be able to survive in the long-term in areas that are not forested, we excluded all individuals 
occurring in cells that were predicted to lose forest cover by 2020 and 2050, respectively. 


QUANTIFICATION AND STATISTICAL ANALYSIS 


As an analytical approach, we used a combination of negative binomial regression models [66] and design-based inference [15, 67] to 
estimate the parameters necessary for building a spatial density distribution model for Bornean orangutans as proposed by Hedley 
et al. [68]. 


Calculating model offset 

In the predictive density distribution model, we used an offset term [69] to convert the number of orangutan nests per transect, into 
the number of individuals per square kilometer. It included the product of the area that was effectively sampled and the relationship 
between number of nests and number of orangutans. The area that was sampled is described by the length of each transect (I) multi- 
plied by twice the ESW. 

The number of orangutans per observed nest was estimated using the proportion of nest builders in a population (p), the daily pro- 
duction rate of nests (r), and the nest decay rate (t), which represents the number of days for which a nest remains visible in the forest 
[13, 14]. For these parameters we used p = 0.88 and r = 1.12 nests/day/individual from Spehar et al. [70], representing a combination 
of the most current nest life-history parameters for Bornean orangutan populations (see below how t was determined). 


Estimation of effective strip width 

For the ground transects, the effective strip width (ESW) was estimated using Distance 6.0 [62]. We used a truncation distance of 
27 m. The models were fitted to the observed data with and without grouping for different habitat categories, using various key 
functions and adjustment terms. The model fit was tested with x? statistics for which we set distance intervals under the “diagnos- 
tics” tab. The fit of the model using habitat specific detection functions was not better than the fit of the model that used a single 
detection function across habitats, as established by Akaike Information Criterion (AIC). As a consequence, we applied a global 
detection function and resulting effective strip width (ESW) to all ground transects. The model with the best fit, based on the lowest 
AIC and x2 statistics, was one with a half-normal key function and a simple polynomial adjustment of order 4. 

Nests with a ppd larger than the truncation distance were excluded from the dataset. We assumed that nests without ppd were 
distributed at similar distances along transects as the nests for which ppds were reported. Therefore, we truncated them by randomly 
excluding the same proportion of nests that were excluded from transects with known distances, leaving 34,415 nests in the dataset. 
The estimated ESW was 15.95 m, and nest detection probabilities for ground transects was 0.59. This is in line with reported detec- 
tion probability for other ape surveys [71]. 

Helicopter surveys did not contain information about the ppds from the transects to the nests. Thus, the ESW for those surveys was 
set to 75 m, which corresponds to half of the maximum visibility from the helicopter to the sides of the survey line [72]. Yearly abun- 
dance estimates were tested for sensitivity to the assumed aerial ESW, but did not vary significantly (abundance estimate with aerial 
ESW = 100 m in comparison to 75 m: correlation coefficient > 0.99, maximum difference 2.127%, aerial ESW = 50 m in comparison to 
75 m: correlation coefficient = 1, maximum difference 3.90496, n = 16 years). 


Estimation of nest decay rate and extrapolation 

We updated the nest decay rate for two sites in the Bornean orangutan range (Sabangau in Central Kalimantan and Lesan in East 
Kalimantan), using the modification of the approach from Laing et al. [66], used in Wich et al. [71]. Additionally, we used site-specific 
decay rates available from the literature for Kinabatangan, Sabah [25] and Gunung Palung, West Kalimantan [63]. For the calculation 
of the nest decay time we used logistic models (left-truncated with normalized intercept, log-transformed and reciprocal) [66] and 
nest age as the only predictor. The product of the daily decay probability and time since nest construction was summed over 
2000 days to calculate mean decay time. The model estimates from the three approaches were model-averaged using their AIC 
weights. The time until nest decay for Sabangau was found to be 496.3 days (n = 423, 95% Cl: 453.1 to 542.9 days) and 
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582.5 days (n = 88, 95% Cl: 461.2 to 753.1) for Lesan, which is similar to the nest decay rate estimated in Spehar et al. [70] for this 
area. We bootstrapped the data 1,000 times and determined the 95% confidence interval by model-averaging the 2.5% and 97.5% 
lower and upper confidence limits. 

The sites, for which we had nest decay values, experience different environmental conditions. The respective values were thus 
used for different parts of the Bornean orangutan range, based on the location of transects within provinces and forest types: (a) 
Sabangau nest decay, 496.3 days (this publication), for peatswamp forests in Central Kalimantan; (b) Lesan nest decay, 583 days 
(this publication), for East and South Kalimantan; (c) Average of Gunung Palung lowland forest, lowland hill and mid-elevation 
nest decay, 276 days [63], for lowland forests in Sarawak, West and Central Kalimantan; (d) Gunung Palung montane forest nest 
decay, 321.3 days [63], for montane forests (> 800 m above sea level (asl)) in Sarawak, West and Central Kalimantan; (e) Gunung 
Palung peatswamp forest nest decay, 399 days [63], for peatswamp forests in West Kalimantan and Sarawak; (f) Kinabatangan 
nest decay, 202 days [25], for Sabah. 


Model structure and multi-model inference 

We used a Generalized Linear Model with a negative binomial error structure and log link function [69] to assess the effect of climate, 
habitat and anthropogenic pressures on orangutans and predict the density distribution across the range. The full model, including all 
predictor variables and the offset term, had the following structure: orangutan nest count on transect ~year + mean temperature + 
rainfall variability + rainfall in dry months + rainfall in dry months? + topsoil organic carbon content + peatswamp cover + lowland 
forest cover + lower montane forest cover + deforestation + human population density + orangutan killing estimates + percent 
population with religious hunting taboos + offset + dispersion parameter. It had been shown that higher orangutan densities occur 
in areas of intermediate levels of rainfall in dry months [11], therefore we included the squared rainfall in dry months. A negative 
coefficient indicates highest orangutan densities at intermediate values of rainfall. 

We tested for collinearity, which was not an issue (largest Variance Inflation Factor = 4.429, see also Table S2) and leverage values 
as well as DFBeta values did not indicate obviously influential cases [73, 74]. The model was not strongly overdispersed (dispersion 
parameter: 1.675). 

As a test of the significance of the predictors, we compared the fit of the full model, as described above, to the null model, only 
including the intercept and the offset term [75]. The comparison was based on a likelihood ratio test. We fitted the models in R 
(version 3.x [56]) using the function glm.nb of the R package MASS and determined Variance Inflation Factors using the function 
vif of the R package car [58]. 

To minimize model uncertainty in spatial model prediction, we applied multi-model inference and assessed all possible combina- 
tions of covariates included in the full model (n = 6,144) [see also 71]. Out of all possible models, only 18 models were in the confi- 
dence set, combining 95% of the AIC weight (Table S3). The best model was the full-model lacking the orangutan killing estimates 
and percent population with religious hunting taboos (Tables S1 and S3). Predictions of all models were averaged, after weighting by 
the models’ AIC weight [76] and used to predict the orangutan density for all 1x1 km cells across the range. We model averaged in link 
space and only after that exponentiated the averaged predictions to get the abundance estimate per grid cell. 

In the output of the density distribution models, all pixels outside the previously defined metapopulations were excluded to avoid 
overestimating Bornean orangutan density, assuming that all larger populations are known to date. Density estimates were summed 
for each metapopulation and land-use category of interest to retrieve total abundance per metapopulation or category [16]. 


Parametric bootstrapping to estimate confidence limits 
The 95% confidence limits of the model predictions were estimated using parametric bootstrapping (n = 1,000). The model-averaged 
fitted estimates and their standard errors (SE), as well as estimate and SE for the dispersion parameter, theta, were used to generate 
1,000 new instances of model estimates by sampling from normal distributions with means and standard deviations being the model 
estimates and their standard errors, respectively. These bootstrapped estimates were then used, together with the model offset and 
the predictors, to sample an instance of the response from a negative binomial distribution with a mean and dispersion parameter 
determined by the bootstrapped estimates. 

We fit the models with the bootstrapped response, resulting in bootstrapped model estimates and AlC-values for each model. 
Using the bootstrapped model-estimates, a prediction was made for each grid cell and study year and from these, the confidence 
limits of the mean and total abundance of cells or groups of cells were determined using the percentile method [77]. 


Spatial overlap of orangutan density distribution and resource use 

With the aim of assessing the differences in the orangutan abundance and change in response to resource use during the survey 
period, we compared the orangutan density distribution from the first and last year of the survey period with maps for land-cover 
classes and area converted into industrial agriculture (oil palm and paper pulp plantations) [32, 55]. The lack of repeat sampling 
through time in areas of land-cover change made it necessary to approach this study in two steps. First, we fitted the model using 
habitat cover and threat predictors and second, overlaid the estimated densities with independent maps of land-cover change to 
infer about patterns of orangutan loss. However, as these maps represent related information, we cannot entirely exclude potential 
circularity in the approach taken. The only approach that completely allows to avoid this problem is to systematically sample across 
gradients of land-use change through time. 
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From the land-use layers we extracted three classes representing changes of orangutan habitat due to resource use (establish- 
ment of industrial oil palm and paper pulp plantations, deforestation, and selective logging) that occurred during the study 
period (1999 — 2015), three classes representing forested areas in 2015 (regrowth forest, primary forest, and primary montane forests 
(> 750 m asl)), and two classes depicting non-forested areas in 2015 (industrial plantations established before 2000 and ‘other’). 
Regrowth forests were areas that were non-forest in 1973, but had forest cover in 2015. The category ‘other’ included scrublands, 
urban, agricultural and non-forest areas that were not contained in the other categories. It was possible that during the study period 
an area was first selectively logged or deforested, and then industrial plantations were established. In our analysis, we counted these 
areas only as industrial plantations, as this was the final stage of the land-use transition. We then pooled the average abundance and 
density in each land-use class or resource use category and calculated the 95% confidence interval. 


DATA AND SOFTWARE AVAILABILITY 


The original orangutan survey, absence and nest decay data used in this study can be requested from the IUCN SSC. A.P.E.S. 
database (http://apesportal.eva.mpg.de/database/archiveTable, ID: MYS IDN Multiple sites 1999 01 01 Voigt Wich et al.). 
Compiled datasets necessary to reproduce the analysis and figures can be downloaded using the link in the database entry. All 
code used for the described analysis is available at https://github.com/MariaVoigt/OU-density-distribution-pipeline.git. 
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Figure S1: Model-averaged Parameter Estimates, Related to Figure 1 and 2. 

Predictors were standardized to a mean of zero and a standard deviation of one to facilitate the comparison. Parameter 
estimates (diamonds) and 95% confidence intervals (horizontal lines) were model averaged, using the AIC weights. The 
interpretation of the linear term of rainfall in dry months depends on whether the quadratic term is included in the model. 


All models which only included the linear (but not the quadratic) term had AIC weights < 0.001. Therefore the coefficient is 
notshown here. 
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Figure S2: Yearly Abundance of All Viable Orangutan Populations between 1999 and 2015, Related to Figure 1. 
Orangutan abundance ofall populations with more than 100 individuals in 2015 over the study period from 1999 to 2015 
including the 95% confidence intervals (grey). Study years are represented on the x-axes. Map shows the location of the 
populations. Abbreviations in population names as follows: Ls = Landscape, frag S = fragmented South. 
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Figure S3: Bornean Orangutan Density and Decline in Resource Use Categories, Related to Figure 3. 

Average density and its change between the first study year (1999) andlast study year (2015) in areas in which industrial oil 
palm or paper pulp plantations were established, which were deforested or selectively logged during the study period (1999 
— 2015), in areas with forest (regrowth, primary and montane forest) and in areas without forest (areas that were transformed 
in plantations before 2000 and otherareas). The “*” indicates the absence of orangutans in the respective category. The error 
bars indicate the 95% confidence interval. The percent orangutan density loss in comparison to 1999 is given in rectangles. 


Table S1: Model Coefficients from Full and Best Model and Summed AIC Weights, Related to Figure 1 and 2. 


Model estimates 
Model AIC weights 


Full Best 
Intercept - -0.56 -0.57 
Year 1.00 -0.22 -0.19 
Mean temperature 1.00 0.80 0.66 
Rainfall variability 1.00 -0.47 -0.49 
Rainfall in dry months 1.00 -1.30 -1.31 
(Rainfall in dry months)? 1.00 -1.08 -1.06 
Topsoil organic carbon content 0.48 -0.08 -0.06 
Peatswamp cover 1.00 0.35 0.35 
Lowland forest cover 1.00 0.52 0.47 
Lower montane forest cover 0.82 -0.14 -0.27 
Deforestation 0.85 0.10 0.09 
Human population density 1.00 -0.53 -0.55 
Orangutan killing estimate 0.14 0.02 - 
Hunting taboo 0.29 0.01 - 
AIC 16204.29 16195.93 
Model weight 0.0037 0.2455 
Model rank (of 6144 models) 25 1 


AIC weights of coefficients are calculated by summing the AIC weights of the models in which the coefficient is present. A 
weight close to 1 indicates an influential predictor. The interpretation of the linear term of rainfall in dry months depends on 
whether the quadratic term is in the model and should not be averaged over all models. All models which only included the 
linear term had AIC weights < 0.001. Their influence on the average coefficient value was thus negligible. For the full and 
the best model the AIC, model weight and model rank are given at the bottom of the table. 


Table S2: Correlation Matrix for the Predictors Used in the Density Distribution Model, Related to Figure 1 and 2. 


Model predictors 





(1) Year 


(2) Mean temperature 


(3) Rainfall variability 


(4) Rainfall in dry months 


(5) Topsoil organic carbon 
content 


(6) Peatswamp cover 


(7) Lowland forest cover 


(8) Lower montane forest 
cover 


(9) Deforestation 


(10) Human population 
density 

(11) Orangutan killing 
estimate 


(12) Hunting taboo 


(1) 


1.000 


0.026 


0.143 


-0.093 


0.066 


0.153 


-0.194 


0.409 


-0.013 


0.052 


0.026 


Q) 


1.000 


0.459 


-0.406 


0.241 


0.219 


0.005 


-0.817 


0.265 


0.338 


-0.51 


0.261 


(3) 


-0.675 


0.223 


0.188 


-0.291 


-0.32 


0.209 


0.484 


-0.254 


0.327 


(4) 


1.000 


-0.118 


-0.129 


0.239 


0.343 


-0.241 


-0.389 


0.289 


-0.602 


(5) 


1.000 


0.370 


-0.204 


-0.116 


0.063 


0.161 


-0.098 


0.022 


(6) 


1.000 


-0.256 


-0.109 


-0.002 


0.063 


-0.043 


0.029 


(7) 


1.000 


-0.296 


-0.362 


-0.468 


0.075 


0.057 


(8) 


1.000 


-0.164 


-0.148 


-0.285 


©) 


1.000 


0.120 


-0.188 


0.122 


(10) 


1.000 


-0.306 


-0.041 


(11) 


1.000 


-0.097 


(12) 


1.000 


Table S3: Models Included in the 95% Confidence Set, Related to Figure 1 and 2. 


AIC 
Model Model coefficients df AIC  AAIC weight 
Nr. Y T RV RD RD? OC PC LC MC DF PD KE HT 
1 14 16195.9 0 045 
2 13 161962 0.3 0.209 
3 15 16198 21 0.086 
4 14 16198.1 22 0.081 
5 14  16198.9 3 0.056 
6 12 16199 3.1 0.054 
7 13 16199.4 3.5 0.044 
8 12 162004 4.5 0.026 
9 13 162004 4.5 0.026 
10 14 16201 5.1 0.019 
11 15 16201 5.1 0.019 
12 13 16201] 52 0.019 
13 14 16201] 52 0.018 
14 13 16202 6.1 0.012 
15 13 162023 6.4 0.01 
16 12 162024 6.5 0.01 
17 13 162024 6.5 0.01 
18 15 162028 6.9 0.008 


The model rank, the included coefficients, degrees of freedom (df), their AIC, AAIC and AIC weights are given. The full 
model included orangutan nest count on transect ~ year (Y) + mean temperature (T) + rainfall variability (RV) + rainfall in 
dry months (RD) + rainfall in dry months? (RD?) + topsoil organic carbon content (OC) + peatswamp cover (PC) + lowland 
forest cover (LC) + lower montane forest cover (MC) + deforestation (DF) + human population density (PD) + orangutan 
killing estimate (KE) + hunting taboo (HT). 
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